/********************************
 * Rachel Robey, Los Alamos High School
 * Copyright 2009
 ********************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

/* Memory allocation routines */
double *dvector(int n);
double **dmatrix(int m, int n);
int **imatrix(int m, int n);

int main(int argc, char *argv[]) {
   double **Values;
   int i, j;
   int h, numpairs;
   double average, vario;
   FILE *fdata, *fsemi;
   char string[80];
   int matrix_size_x, matrix_size_y;
   matrix_size_x = 50;
   matrix_size_y = 50;

   Values = dmatrix(matrix_size_y, matrix_size_x);

   //Read in Data
   fdata = fopen("sgsim.out", "r");
   for(j = 0; j < matrix_size_y; j++) {
      for(i = 0; i < matrix_size_x; i++) {
         fgets(string, 80, fdata);
         sscanf(string, "%lf", &Values[j][i]);
      }
   }
   fsemi = fopen("data.csv", "w");
//Semivariogram
   //X Direction
   for(h = 0; h <= matrix_size_x-1; h++) {
      numpairs = 0;
      average = 0;
      for(j = 0; j < matrix_size_y; j++) {
         for(i = 0; i < matrix_size_x-h; i++) {
            if(Values[j][i] >= 0 && Values[j][i+h] >= 0) {
               average += (Values[j][i]-Values[j][i+h])*(Values[j][i]-Values[j][i+h]);
               numpairs += 1;
            }
         }
      }
      vario = average/(2*numpairs);
      fprintf(fsemi, "%d, %lf\n", h, vario);
   }
   //y Direction
   for(h = 0; h <= matrix_size_y-1; h++) {
      numpairs = 0;
      average = 0;
      for(i = 0; i < matrix_size_x; i++) {
         for(j = 0; j < matrix_size_y-h; j++) {
            if(Values[j][i] >= 0 && Values[j+h][i] >= 0) {
               average += (Values[j][i]-Values[j+h][i])*(Values[j][i]-Values[j+h][i]);
               numpairs += 1;
            }
         }
      }
      vario = average/(2*numpairs);
      fprintf(fsemi, "%d, %lf\n", h, vario);
   }
   exit(0);
}

